Rice-Wheat System Spatial Decision Support Tool

1 Introduction

We present a rice-wheat systen spatial decisionb support tool which relies on data from the matched landscape crop assessment surveys in eastern India and a multivariate geoadditive Bayesian model to predict entry points for system optimization in the area of interest.

2 Crop specific regressions

2.1 Rice

library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# library(mvnchol)
library(BayesX)
library(R2BayesX)
library(sf)
library(spdep)
library(rio)
Irrig_Rev_rice_wheat <- import("data/Irrig_Rev_rice_wheat.csv")

library(janitor)
# Irrig_Rev_rice_wheat <- clean_names(Irrig_Rev_rice_wheat)
library(tidyr)
#Irrig_Rev_rice_wheat <- Irrig_Rev_rice_wheat %>% drop_na()

# library(lubridate)
# library(anytime)
# Irrig_Rev_rice_wheat$januaryfirst2017 <- ymd("2017-01-01")

# Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day <- Irrig_Rev_rice_wheat$sowdate_fmt - Irrig_Rev_rice_wheat$januaryfirst2017

# Irrig_Rev_rice_wheat$sowdate_fmt_rice_day <- Irrig_Rev_rice_wheat$sowdate_fmt_rice - Irrig_Rev_rice_wheat$januaryfirst2017

Irrig_Rev_rice_wheat$harvest_day_rice <- Irrig_Rev_rice_wheat$l_crop_duration_days_rice + Irrig_Rev_rice_wheat$sowdate_fmt_rice_day

# Irrig_Rev_rice_wheat$harvest_day_wheat <- Irrig_Rev_rice_wheat$l_crop_duration_days + Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day


shpname <- file.path(getwd(), "shp", "India_aoi_sf_sp")
India_aoi_sp_bnd <- BayesX::shp2bnd(shpname = shpname, regionnames = "District", check.is.in = F)

Reading map ... finished
Note: map consists originally of 50 polygons
Note: map consists of 47 regions

f_rice_yield_MRF <- list(
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)



K <- neighbormatrix(India_aoi_sp_bnd)
head(K)

           Araria Arwal Aurangabad Banka Begusarai Bhagalpur Arah Buxar
Araria          4     0          0     0         0         0    0     0
Arwal           0     6         -1     0         0         0   -1     0
Aurangabad      0    -1          3     0         0         0    0     0
Banka           0     0          0     3         0        -1    0     0
Begusarai       0     0          0     0         5         0    0     0
Bhagalpur       0     0          0    -1         0         6    0     0
           Darbhanga Gaya Gopalganj Jamui Jehanabad Kaimur Katihar Khagaria
Araria             0    0         0     0         0      0       0        0
Arwal              0   -1         0     0        -1      0       0        0
Aurangabad         0   -1         0     0         0      0       0        0
Banka              0    0         0    -1         0      0       0        0
Begusarai          0    0         0     0         0      0       0       -1
Bhagalpur          0    0         0     0         0      0      -1       -1
           Kishanganj Lakhisarai Madhepura Madhubani Munger Muzaffarpur Nalanda
Araria             -1          0        -1         0      0           0       0
Arwal               0          0         0         0      0           0       0
Aurangabad          0          0         0         0      0           0       0
Banka               0          0         0         0     -1           0       0
Begusarai           0         -1         0         0     -1           0       0
Bhagalpur           0          0        -1         0     -1           0       0
           Nawada WestChamparan Patna EastChamparan Purnia Rohtas Saharsa
Araria          0             0     0             0     -1      0       0
Arwal           0             0    -1             0      0     -1       0
Aurangabad      0             0     0             0      0     -1       0
Banka           0             0     0             0      0      0       0
Begusarai       0             0    -1             0      0      0       0
Bhagalpur       0             0     0             0     -1      0       0
           Samastipur Saran Sheikhpura Sheohar Sitamarhi Siwan Supaul Vaishali
Araria              0     0          0       0         0     0     -1        0
Arwal               0     0          0       0         0     0      0        0
Aurangabad          0     0          0       0         0     0      0        0
Banka               0     0          0       0         0     0      0        0
Begusarai          -1     0          0       0         0     0      0        0
Bhagalpur           0     0          0       0         0     0      0        0
           Balia Chandauli Deoria Gazipur Gorakhpur Kushinagar Maharajganj Mau
Araria         0         0      0       0         0          0           0   0
Arwal          0         0      0       0         0          0           0   0
Aurangabad     0         0      0       0         0          0           0   0
Banka          0         0      0       0         0          0           0   0
Begusarai      0         0      0       0         0          0           0   0
Bhagalpur      0         0      0       0         0          0           0   0
           Siddharthnagar
Araria                  0
Arwal                   0
Aurangabad              0
Banka                   0
Begusarai               0
Bhagalpur               0

## Also need to transform to factor for
## setting up the MRF smooth.
Irrig_Rev_rice_wheat$District <- as.factor(Irrig_Rev_rice_wheat$a_q103_district)

## Now note that not all regions are observed,
## therefore we need to remove those regions
## from the penalty matrix
rn <- rownames(K)
lv <- levels(Irrig_Rev_rice_wheat$District)
i <- rn %in% lv
K <- K[i, i]

set.seed(321)
library(bamlss)
b_rice_yield_MRF <- bamlss(f_rice_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")

AICc 12433.92 logPost -6216.75 logLik -6063.08 edf 148.79 eps 2.3399 iteration   1
AICc 12239.92 logPost -5910.29 logLik -5967.69 edf 147.28 eps 0.4613 iteration   2
AICc 12211.69 logPost -5884.01 logLik -5956.93 edf 144.14 eps 0.1461 iteration   3
AICc 12205.11 logPost -5882.70 logLik -5953.35 edf 144.42 eps 2.9473 iteration   4
AICc 12203.12 logPost -5882.46 logLik -5952.28 edf 144.48 eps 0.0429 iteration   5
AICc 12202.50 logPost -5882.43 logLik -5951.95 edf 144.50 eps 0.0113 iteration   6
AICc 12202.32 logPost -5882.45 logLik -5951.84 edf 144.52 eps 0.0088 iteration   7
AICc 12202.28 logPost -5882.44 logLik -5951.82 edf 144.52 eps 0.0010 iteration   8
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0002 iteration   9
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration  10
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration  11
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration  12
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration  12
elapsed time:  9.60sec
Starting the sampler...

|                    |   0%  1.84min
|*                   |   5%  1.80min  5.70sec
|**                  |  10%  1.74min 11.60sec
|***                 |  15%  1.61min 17.10sec
|****                |  20%  1.54min 23.10sec
|*****               |  25%  1.48min 29.53sec
|******              |  30%  1.40min 35.96sec
|*******             |  35%  1.30min 41.90sec
|********            |  40%  1.18min 47.01sec
|*********           |  45%  1.06min 51.85sec
|**********          |  50% 56.70sec 56.70sec
|***********         |  55% 50.38sec  1.03min
|************        |  60% 44.43sec  1.11min
|*************       |  65% 38.66sec  1.20min
|**************      |  70% 32.88sec  1.28min
|***************     |  75% 27.22sec  1.36min
|****************    |  80% 21.67sec  1.44min
|*****************   |  85% 16.16sec  1.53min
|******************  |  90% 10.72sec  1.61min
|******************* |  95%  5.34sec  1.69min
|********************| 100%  0.00sec  1.78min

## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_rice_yield_MRF)

Call:
bamlss(formula = f_rice_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian 
Link function: mu = identity, sigma = log
*---
Formula mu:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.90903 0.53256 1.64495 3.75492      0.013
rice_duration_class_long 0.07427 0.01184 0.07438 0.13788      0.075
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(sowdate_fmt_rice_day).tau21     5.510e-01 1.460e-02 2.985e-01 2.446e+00
s(sowdate_fmt_rice_day).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf       3.784e+00 1.625e+00 3.700e+00 5.909e+00
s(g_q5305_irrig_times_rice).tau21 2.324e+00 2.653e-01 1.653e+00 8.760e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.482e+00 4.641e+00 6.530e+00 7.994e+00
s(nperha_rice).tau21              4.816e-01 1.694e-04 7.272e-02 3.602e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                2.866e+00 1.010e+00 2.509e+00 6.675e+00
s(p2o5perha_rice).tau21           9.498e-02 9.163e-05 7.279e-03 8.451e-01
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             1.889e+00 1.007e+00 1.411e+00 4.674e+00
s(District,id='mrf1').tau21       4.290e-02 7.425e-05 1.089e-02 2.036e-01
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         2.274e+01 1.975e+00 2.845e+01 3.410e+01
s(District,id='re2').tau21        5.664e+00 2.034e-01 5.310e+00 1.506e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.560e+01 3.391e+01 3.581e+01 3.594e+01
                                  parameters
s(sowdate_fmt_rice_day).tau21          0.182
s(sowdate_fmt_rice_day).alpha             NA
s(sowdate_fmt_rice_day).edf            3.291
s(g_q5305_irrig_times_rice).tau21      2.895
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        7.106
s(nperha_rice).tau21                   1.398
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.489
s(p2o5perha_rice).tau21                0.018
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  1.787
s(District,id='mrf1').tau21            0.029
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             31.761
s(District,id='re2').tau21            12.212
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.917
---
Formula sigma:
---
sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + 
    s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                              Mean      2.5%       50%     97.5% parameters
(Intercept)              -0.119970 -0.179496 -0.119649 -0.063475     -0.114
rice_duration_class_long  0.061685  0.009405  0.061474  0.115939      0.060
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9860 0.8949 1.0000     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(sowdate_fmt_rice_day).tau21     2.255e-02 9.588e-05 4.358e-03 1.617e-01
s(sowdate_fmt_rice_day).alpha     9.779e-01 8.326e-01 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf       1.623e+00 1.009e+00 1.340e+00 3.522e+00
s(g_q5305_irrig_times_rice).tau21 2.200e-01 2.578e-04 1.086e-01 1.030e+00
s(g_q5305_irrig_times_rice).alpha 9.450e-01 6.660e-01 9.945e-01 1.000e+00
s(g_q5305_irrig_times_rice).edf   4.016e+00 1.123e+00 4.209e+00 6.532e+00
s(nperha_rice).tau21              1.741e-01 4.189e-04 3.295e-02 1.329e+00
s(nperha_rice).alpha              9.529e-01 6.319e-01 9.979e-01 1.000e+00
s(nperha_rice).edf                2.684e+00 1.043e+00 2.352e+00 6.037e+00
s(p2o5perha_rice).tau21           3.062e-01 1.414e-04 1.220e-02 2.836e+00
s(p2o5perha_rice).alpha           9.579e-01 6.586e-01 9.982e-01 1.000e+00
s(p2o5perha_rice).edf             2.461e+00 1.017e+00 1.808e+00 6.536e+00
s(District,id='mrf1').tau21       2.374e-03 7.354e-05 1.646e-03 7.531e-03
s(District,id='mrf1').alpha       8.160e-01 2.828e-01 9.140e-01 1.000e+00
s(District,id='mrf1').edf         1.912e+01 2.882e+00 2.035e+01 2.921e+01
s(District,id='re2').tau21        2.233e-02 5.515e-03 2.132e-02 4.458e-02
s(District,id='re2').alpha        7.382e-01 1.551e-01 8.258e-01 1.000e+00
s(District,id='re2').edf          2.842e+01 1.981e+01 2.919e+01 3.199e+01
                                  parameters
s(sowdate_fmt_rice_day).tau21          0.014
s(sowdate_fmt_rice_day).alpha             NA
s(sowdate_fmt_rice_day).edf            1.802
s(g_q5305_irrig_times_rice).tau21      0.191
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        4.768
s(nperha_rice).tau21                   0.091
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     3.139
s(p2o5perha_rice).tau21                2.253
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  6.288
s(District,id='mrf1').tau21            0.007
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             28.648
s(District,id='re2').tau21             0.002
s(District,id='re2').alpha                NA
s(District,id='re2').edf              10.540
---
Sampler summary:
-
DIC = 12123.05 logLik = -6011.865 pd = 99.3237
runtime = 107.63
---
Optimizer summary:
-
AICc = 12202.27 edf = 144.5353 logLik = -5951.809
logPost = -5882.444 nobs = 4537 runtime = 9.6

# Plot the nonlinear effect
plot(b_rice_yield_MRF, model = "mu", term = "s(sowdate_fmt_rice_day)")

plot(b_rice_yield_MRF, model = "mu", term = "s(g_q5305_irrig_times_rice)")

plot(b_rice_yield_MRF, model = "mu", term = "s(nperha_rice)")

plot(b_rice_yield_MRF, model = "mu", term = "s(p2o5perha_rice)")

2.2 Wheat model

f_wheat_yield_MRF <- list(
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)
set.seed(321)
b_wheat_yield_MRF <- bamlss(f_wheat_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")

AICc 8958.541 logPost -4419.70 logLik -4335.18 edf 139.61 eps 0.3246 iteration   1
AICc 8445.144 logPost -4096.84 logLik -4071.93 edf 145.76 eps 0.1431 iteration   2
AICc 8344.698 logPost -4042.08 logLik -4030.62 edf 137.40 eps 0.0345 iteration   3
AICc 8319.455 logPost -4032.32 logLik -4022.32 edf 133.33 eps 0.0132 iteration   4
AICc 8313.001 logPost -4030.95 logLik -4020.60 edf 131.91 eps 0.0050 iteration   5
AICc 8312.046 logPost -4029.90 logLik -4020.23 edf 131.81 eps 0.0019 iteration   6
AICc 8311.784 logPost -4029.54 logLik -4020.15 edf 131.76 eps 0.0010 iteration   7
AICc 8311.335 logPost -4029.64 logLik -4020.12 edf 131.57 eps 0.0005 iteration   8
AICc 8311.297 logPost -4029.53 logLik -4020.11 edf 131.57 eps 0.0003 iteration   9
AICc 8311.268 logPost -4029.43 logLik -4020.10 edf 131.56 eps 0.0002 iteration  10
AICc 8311.245 logPost -4029.35 logLik -4020.09 edf 131.56 eps 0.0002 iteration  11
AICc 8311.226 logPost -4029.28 logLik -4020.09 edf 131.56 eps 0.0002 iteration  12
AICc 8311.211 logPost -4029.22 logLik -4020.08 edf 131.56 eps 0.0001 iteration  13
AICc 8311.197 logPost -4029.16 logLik -4020.08 edf 131.55 eps 0.0001 iteration  14
AICc 8311.185 logPost -4029.11 logLik -4020.07 edf 131.55 eps 0.0001 iteration  15
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration  16
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration  16
elapsed time:  6.08sec
Starting the sampler...

|                    |   0%  1.53min
|*                   |   5%  1.38min  4.35sec
|**                  |  10%  1.32min  8.83sec
|***                 |  15%  1.24min 13.11sec
|****                |  20%  1.18min 17.67sec
|*****               |  25%  1.14min 22.74sec
|******              |  30%  1.05min 26.92sec
|*******             |  35% 57.61sec 31.02sec
|********            |  40% 52.11sec 34.74sec
|*********           |  45% 46.87sec 38.35sec
|**********          |  50% 42.02sec 42.02sec
|***********         |  55% 37.54sec 45.88sec
|************        |  60% 32.99sec 49.49sec
|*************       |  65% 28.57sec 53.05sec
|**************      |  70% 24.65sec 57.52sec
|***************     |  75% 23.26sec  1.16min
|****************    |  80% 21.38sec  1.43min
|*****************   |  85% 16.78sec  1.59min
|******************  |  90% 11.20sec  1.68min
|******************* |  95%  5.54sec  1.75min
|********************| 100%  0.00sec  1.83min

## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_wheat_yield_MRF)

Call:
bamlss(formula = f_wheat_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian 
Link function: mu = identity, sigma = log
*---
Formula mu:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.5740916 -1.1340813 -0.6630343  0.1609182     -1.359
variety_type_NMWV    0.3056829  0.2611116  0.3052558  0.3529448      0.303
g_q5305_irrig_times  0.4153445  0.3895350  0.4153230  0.4439696      0.411
nperha               0.0012336  0.0006668  0.0012288  0.0017957      0.001
p2o5perha            0.0031334  0.0020527  0.0031371  0.0042060      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                    Mean      2.5%       50%     97.5%
s(sowdate_fmt_wheat_day).tau21 4.037e-02 7.694e-05 3.583e-03 3.627e-01
s(sowdate_fmt_wheat_day).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_wheat_day).edf   1.940e+00 1.015e+00 1.491e+00 4.943e+00
s(District,id='mrf1').tau21    3.294e-03 6.895e-05 9.061e-04 1.507e-02
s(District,id='mrf1').alpha    1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf      1.865e+01 3.933e+00 1.853e+01 3.206e+01
s(District,id='re2').tau21     4.793e+00 1.759e+00 4.533e+00 9.645e+00
s(District,id='re2').alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf       3.589e+01 3.577e+01 3.590e+01 3.596e+01
                               parameters
s(sowdate_fmt_wheat_day).tau21     15.401
s(sowdate_fmt_wheat_day).alpha         NA
s(sowdate_fmt_wheat_day).edf        8.309
s(District,id='mrf1').tau21         0.003
s(District,id='mrf1').alpha            NA
s(District,id='mrf1').edf          26.415
s(District,id='re2').tau21          6.085
s(District,id='re2').alpha             NA
s(District,id='re2').edf           35.897
---
Formula sigma:
---
sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + 
    nperha + p2o5perha + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.7819694 -0.8953612 -0.7812289 -0.6696800     -0.316
variety_type_NMWV    0.2056530  0.1502403  0.2059671  0.2547450      0.198
g_q5305_irrig_times  0.0927109  0.0623800  0.0929975  0.1245421      0.098
nperha              -0.0011626 -0.0018926 -0.0011700 -0.0004521     -0.001
p2o5perha            0.0017864  0.0003726  0.0017712  0.0031635      0.002
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.8957 0.4855 0.9745     1
-
Smooth terms:
                                    Mean      2.5%       50%     97.5%
s(sowdate_fmt_wheat_day).tau21 8.472e-02 5.798e-04 4.162e-02 4.208e-01
s(sowdate_fmt_wheat_day).alpha 9.647e-01 7.699e-01 9.984e-01 1.000e+00
s(sowdate_fmt_wheat_day).edf   2.574e+00 1.067e+00 2.524e+00 4.535e+00
s(District,id='mrf1').tau21    6.406e-04 5.097e-05 3.717e-04 2.418e-03
s(District,id='mrf1').alpha    8.975e-01 4.521e-01 9.764e-01 1.000e+00
s(District,id='mrf1').edf      1.089e+01 2.084e+00 9.704e+00 2.304e+01
s(District,id='re2').tau21     1.986e-02 9.452e-03 1.872e-02 3.562e-02
s(District,id='re2').alpha     7.360e-01 1.839e-01 8.280e-01 1.000e+00
s(District,id='re2').edf       2.834e+01 2.429e+01 2.855e+01 3.129e+01
                               parameters
s(sowdate_fmt_wheat_day).tau21      0.035
s(sowdate_fmt_wheat_day).alpha         NA
s(sowdate_fmt_wheat_day).edf        2.415
s(District,id='mrf1').tau21         0.001
s(District,id='mrf1').alpha            NA
s(District,id='mrf1').edf          13.967
s(District,id='re2').tau21          0.203
s(District,id='re2').alpha             NA
s(District,id='re2').edf           34.551
---
Sampler summary:
-
DIC = 8219.792 logLik = -4070.916 pd = 77.9595
runtime = 110.45
---
Optimizer summary:
-
AICc = 8311.175 edf = 131.5535 logLik = -4020.075
logPost = -4029.064 nobs = 4537 runtime = 6.08

# Plot the nonlinear effect
plot(b_wheat_yield_MRF, model = "mu", term = "s(sowdate_fmt_wheat_day)")

## Now, to predict the spatial effects we set up new data.
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))

## Predict for the structured spatial effects.
p_str_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)

## And the unstructured spatial effect.
p_unstr_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)

## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
    x = p_str_wheat_yield_MRF$mu, id = nd$District,
    main = expression(mu), title = "Structured spatial effect"
)

plotmap(India_aoi_sp_bnd,
    x = p_str_wheat_yield_MRF$sigma, id = nd$District,
    main = expression(sigma), title = "Structured spatial effect"
)

## Random effects.
plotmap(India_aoi_sp_bnd,
    x = p_unstr_wheat_yield_MRF$mu, id = nd$District, title = "Unstructured spatial effect"
)

plotmap(India_aoi_sp_bnd,
    x = p_unstr_wheat_yield_MRF$sigma, id = nd$District, title = "Unstructured spatial effect"
)

3 Multivariate Spatial Geoadditive Bayesian Regression Model

4 Spatial + for multivariate model

library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# install_github("https://github.com/meteosimon/mvnchol")
library(mvnchol)

# Remove NAs in the monsoon variables
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$onset_2017)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$monsoon_onset_dev)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$median_onset_82_15)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$sd_onset_82_15)))

# Rice first stage
f_sow_rice_1st_stage <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_rice_1st_stage_MRF <- bamlss(f_sow_rice_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")

AICc 38862.03 logPost -165306. logLik -19371.2 edf 59.010 eps 0.1226 iteration   1
AICc 34841.08 logPost -31072.7 logLik -17339.3 edf 79.788 eps 0.0980 iteration   2
AICc 33498.94 logPost -18399.0 logLik -16648.2 edf 98.970 eps 0.0632 iteration   3
AICc 33180.98 logPost -17093.4 logLik -16491.0 edf 97.287 eps 0.0297 iteration   4
AICc 33096.83 logPost -16986.5 logLik -16452.6 edf 93.751 eps 0.0069 iteration   5
AICc 33056.53 logPost -16935.2 logLik -16435.7 edf 90.674 eps 0.0017 iteration   6
AICc 33040.37 logPost -16899.0 logLik -16429.8 edf 88.587 eps 0.0009 iteration   7
AICc 33035.63 logPost -16869.0 logLik -16428.1 edf 87.868 eps 0.0003 iteration   8
AICc 33028.91 logPost -16845.9 logLik -16427.6 edf 85.110 eps 0.0001 iteration   9
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration  10
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration  10
elapsed time:  3.76sec
Starting the sampler...

|                    |   0%  1.45min
|*                   |   5%  1.03min  3.26sec
|**                  |  10% 58.05sec  6.45sec
|***                 |  15% 54.63sec  9.64sec
|****                |  20% 52.12sec 13.03sec
|*****               |  25% 49.77sec 16.59sec
|******              |  30% 47.53sec 20.37sec
|*******             |  35% 44.31sec 23.86sec
|********            |  40% 40.83sec 27.22sec
|*********           |  45% 37.08sec 30.34sec
|**********          |  50% 33.14sec 33.14sec
|***********         |  55% 29.32sec 35.84sec
|************        |  60% 25.67sec 38.50sec
|*************       |  65% 22.12sec 41.08sec
|**************      |  70% 18.71sec 43.66sec
|***************     |  75% 15.46sec 46.37sec
|****************    |  80% 12.27sec 49.08sec
|*****************   |  85%  9.12sec 51.69sec
|******************  |  90%  6.02sec 54.16sec
|******************* |  95%  2.99sec 56.80sec
|********************| 100%  0.00sec 59.22sec

fitted_f_sow_rice_1st_stage_MRF <- f_sow_rice_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_rice_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_rice_day - fitted_f_sow_rice_1st_stage_MRF$mu

# Wheat first stage
f_sow_wheat_1st_stage <- list(
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(gw_2018) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_wheat_1st_stage_MRF <- bamlss(f_sow_wheat_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")

AICc 49188.65 logPost -488850. logLik -24543.4 edf 50.285 eps 0.2566 iteration   1
AICc 44824.42 logPost -64385.9 logLik -22386.0 edf 26.046 eps 0.0753 iteration   2
AICc 40409.05 logPost -24539.6 logLik -20136.0 edf 67.444 eps 0.0824 iteration   3
AICc 36648.46 logPost -18973.1 logLik -18246.6 edf 76.219 eps 0.0789 iteration   4
AICc 33900.77 logPost -17313.2 logLik -16871.9 edf 77.104 eps 0.0796 iteration   5
AICc 32636.01 logPost -16736.4 logLik -16238.1 edf 78.490 eps 0.0660 iteration   6
AICc 32419.72 logPost -16629.9 logLik -16129.2 edf 79.163 eps 0.0324 iteration   7
AICc 32413.31 logPost -16627.4 logLik -16125.7 edf 79.502 eps 0.0056 iteration   8
AICc 32412.92 logPost -16627.3 logLik -16125.4 edf 79.574 eps 0.0001 iteration   9
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration  10
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration  10
elapsed time:  1.96sec
Starting the sampler...

|                    |   0% 36.89sec
|*                   |   5% 36.29sec  1.91sec
|**                  |  10% 33.75sec  3.75sec
|***                 |  15% 31.28sec  5.52sec
|****                |  20% 30.00sec  7.50sec
|*****               |  25% 28.23sec  9.41sec
|******              |  30% 26.76sec 11.47sec
|*******             |  35% 25.48sec 13.72sec
|********            |  40% 24.33sec 16.22sec
|*********           |  45% 23.47sec 19.20sec
|**********          |  50% 21.86sec 21.86sec
|***********         |  55% 19.92sec 24.35sec
|************        |  60% 17.92sec 26.88sec
|*************       |  65% 15.81sec 29.36sec
|**************      |  70% 13.71sec 32.00sec
|***************     |  75% 11.54sec 34.61sec
|****************    |  80%  9.30sec 37.20sec
|*****************   |  85%  7.01sec 39.74sec
|******************  |  90%  4.72sec 42.49sec
|******************* |  95%  2.37sec 45.00sec
|********************| 100%  0.00sec 47.50sec

fitted_f_sow_wheat_1st_stage_MRF <- f_sow_wheat_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_wheat_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day - fitted_f_sow_wheat_1st_stage_MRF$mu


# Multivariate with sowing dates as endogenous
f_rice_wheat_yield_MRF_corr_endo <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo <- bamlss(f_rice_wheat_yield_MRF_corr_endo, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)

AICc 349476.1 logPost -2543606 logLik -174262. edf 430.14 eps 1.0000 iteration   1
AICc 157743.5 logPost -2447306 logLik -78382.3 edf 441.55 eps 0.1825 iteration   2
AICc 104307.5 logPost -2420173 logLik -51667.3 edf 439.10 eps 0.2070 iteration   3
AICc 89316.74 logPost -2412636 logLik -44189.9 edf 424.43 eps 0.1254 iteration   4
AICc 86452.65 logPost -2411216 logLik -42768.9 edf 415.34 eps 0.0628 iteration   5
AICc 86231.98 logPost -2411101 logLik -42661.9 edf 412.53 eps 0.0202 iteration   6
AICc 86215.73 logPost -2411097 logLik -42659.1 edf 408.13 eps 0.0027 iteration   7
AICc 86213.14 logPost -2411094 logLik -42658.2 edf 407.83 eps 0.0007 iteration   8
AICc 86211.64 logPost -2411092 logLik -42657.7 edf 407.62 eps 0.0003 iteration   9
AICc 86210.65 logPost -2411090 logLik -42657.4 edf 407.43 eps 0.0002 iteration  10
AICc 86205.44 logPost -2411208 logLik -42657.3 edf 405.37 eps 0.0001 iteration  11
AICc 86205.08 logPost -2411207 logLik -42657.3 edf 405.26 eps 0.0001 iteration  12
AICc 86202.45 logPost -2412410 logLik -42657.2 edf 404.20 eps 0.0001 iteration  13
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration  14
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration  14
elapsed time: 52.51sec
Starting the sampler...

|                    |   0%  8.77min
|*                   |   5%  7.86min 24.83sec
|**                  |  10%  9.20min  1.02min
|***                 |  15%  8.52min  1.50min
|****                |  20%  7.77min  1.94min
|*****               |  25%  7.19min  2.40min
|******              |  30%  7.00min  3.00min
|*******             |  35%  6.43min  3.46min
|********            |  40%  6.07min  4.05min
|*********           |  45%  5.54min  4.53min
|**********          |  50%  4.97min  4.97min
|***********         |  55%  4.55min  5.56min
|************        |  60%  4.15min  6.22min
|*************       |  65%  3.63min  6.75min
|**************      |  70%  3.16min  7.38min
|***************     |  75%  2.69min  8.07min
|****************    |  80%  2.19min  8.77min
|*****************   |  85%  1.66min  9.43min
|******************  |  90%  1.13min 10.20min
|******************* |  95% 34.28sec 10.85min
|********************| 100%  0.00sec 11.74min

summary(multivariate_geo_sow_MRF_corr_endo)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              68.6313 65.0608 69.0601 70.4619      2.030
rice_duration_class_long -1.3633 -2.0104 -1.3634 -0.7134     -1.341
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            2.095e+00 1.260e-04 1.210e-02 2.125e+01    428.875
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.326e+00 9.987e-01 1.008e+00 3.388e+00      7.459
s(onset_2017).tau21         7.779e+00 1.169e-04 2.620e-01 7.322e+01   3370.678
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.790e+00 9.986e-01 1.275e+00 4.686e+00      8.245
s(monsoon_onset_dev).tau21  7.665e+00 2.474e-04 1.551e-01 7.089e+01     16.446
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    1.930e+00 9.996e-01 1.214e+00 5.147e+00      3.687
s(median_onset_82_15).tau21 2.959e+01 8.439e-04 1.024e+01 1.830e+02     18.779
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   3.182e+00 1.001e+00 3.276e+00 6.208e+00      3.766
s(sd_onset_82_15).tau21     1.316e+01 1.644e-03 1.563e+00 8.969e+01      0.537
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       2.170e+00 1.002e+00 1.801e+00 5.006e+00      1.379
s(District,id='mrf1').tau21 3.089e+01 1.575e+01 2.914e+01 5.586e+01   1484.368
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.430e+01 3.394e+01 3.431e+01 3.457e+01     34.979
s(District,id='re2').tau21  1.606e+04 1.019e+04 1.539e+04 2.517e+04      1.087
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.599e+01 3.599e+01 3.600e+01     35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              2.16281 0.43171 1.97757 3.90909      0.000
rice_duration_class_long 0.10345 0.03601 0.10329 0.17645      0.096
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             2.939e+00 7.670e-02 1.377e+00 1.492e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.543e+00 2.630e+00 5.592e+00 8.087e+00
s(g_q5305_irrig_times_rice).tau21 1.715e+00 2.310e-01 1.227e+00 5.434e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.171e+00 4.441e+00 6.195e+00 7.609e+00
s(nperha_rice).tau21              7.856e-01 5.357e-03 3.042e-01 3.674e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                3.921e+00 1.288e+00 3.786e+00 6.668e+00
s(p2o5perha_rice).tau21           2.061e-01 1.334e-04 3.188e-02 1.348e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             2.348e+00 1.010e+00 2.059e+00 5.205e+00
s(District,id='mrf1').tau21       3.043e-03 8.435e-05 1.067e-03 1.628e-02
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         1.462e+01 2.134e+00 1.342e+01 3.017e+01
s(District,id='re2').tau21        4.974e+00 1.229e-01 4.138e+00 1.460e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.529e+01 3.323e+01 3.577e+01 3.593e+01
                                  parameters
s(Res_rice_sow).tau21                 16.244
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.170
s(g_q5305_irrig_times_rice).tau21      1.698
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.538
s(nperha_rice).tau21                   0.991
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.103
s(p2o5perha_rice).tau21                0.172
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.239
s(District,id='mrf1').tau21            0.035
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             32.259
s(District,id='re2').tau21             8.948
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.894
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       114.655 109.483 114.481 120.587      3.830
variety_type_NMWV  -3.072  -3.742  -3.074  -2.460     -3.073
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     228.435    44.388   160.237   704.149   1035.210
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.934     4.524     5.898     7.396      8.391
s(District,id='mrf1').tau21    54.003     9.448    45.568   154.095   4326.923
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.590    33.755    34.685    34.847     34.994
s(District,id='re2').tau21  49904.116 32148.970 48241.930 80332.511      1.087
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.999    35.998    35.999    35.999     35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.3451667 -0.9637935 -0.4133652  0.3762429     -1.419
variety_type_NMWV    0.3422512  0.2824899  0.3429785  0.3966410      0.337
g_q5305_irrig_times  0.4219196  0.3932712  0.4221273  0.4494278      0.424
nperha               0.0015496  0.0009547  0.0015550  0.0021134      0.002
p2o5perha            0.0025587  0.0014420  0.0025588  0.0036574      0.002
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      1.861e-02 7.837e-05 1.649e-03 1.452e-01      0.238
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.607e+00 1.011e+00 1.202e+00 3.875e+00      4.795
s(District,id='mrf1').tau21 1.002e-02 5.295e-05 5.635e-03 4.258e-02      0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   2.528e+01 2.963e+00 2.938e+01 3.356e+01     31.282
s(District,id='re2').tau21  3.750e+00 1.161e+00 3.640e+00 7.588e+00      5.082
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.586e+01 3.568e+01 3.588e+01 3.594e+01     35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.216 -2.236 -2.216 -2.194     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9896 0.9035 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06651 0.04523 0.06601 0.08788      0.074
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9913 0.9196 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.145 -2.167 -2.145 -2.126     -2.141
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9896 0.9163 1.0000     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5178 0.4983 0.5173 0.5394      0.527
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9910 0.9175 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept) -0.011213 -0.045104 -0.007393  0.028365     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.007309 0.003664 0.007292 0.010696      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.01113 0.00202 0.01088 0.01917      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023714 -0.007921  0.023630  0.055818      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) -0.2382 -0.2908 -0.2359 -0.1996     -0.036
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(District,id='mrf1').tau21 2.387e-03 1.004e-04 2.022e-03 6.647e-03      0.005
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.479e+01 1.863e+00 1.601e+01 2.430e+01     22.153
s(District,id='re2').tau21  9.553e-03 1.586e-04 6.116e-03 3.340e-02      0.024
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    1.374e+01 6.443e-01 1.415e+01 2.686e+01     24.454
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.07590 -0.11822 -0.08111 -0.01156     -0.009
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 85815.46 logLik = -42797.02 pd = 221.4134
runtime = 706.63
---
Optimizer summary:
-
AICc = 86202.37 edf = 404.1805 logLik = -42657.27
logPost = -2412333 nobs = 4527 runtime = 52.51

nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))

# Focusing on the cross-equation correlations

## Predict for the structured spatial effects.
p_str_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)

p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt <- as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y)


## And the unstructured spatial effect.
p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)

# Rice sowing spatial equation

plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Structured spatial effect")

plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Unstructured spatial effect")

# Rice yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Structured spatial effect")

plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Unstructured spatial effect")

# Wheat sowing equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Structured spatial effect")

plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Unstructured spatial effect")

# Wheat yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Structured spatial effect")

plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Unstructured spatial effect")

# Focusing on the cross-equation correlations

## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
    x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District,
    main = expression(lambda(rice, wheat)), title = "Structured spatial effect"
)

## Random effects.
plotmap(India_aoi_sp_bnd,
    x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District, main = expression(lambda(rice, wheat)), title = "Unstructured spatial effect"
)

# Rice sowing equation : Non-linear relationships
# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(gw_2017)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(onset_2017)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(monsoon_onset_dev)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(median_onset_82_15)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(sd_onset_82_15)")

# Rice yield equation
# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(Res_rice_sow)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(g_q5305_irrig_times_rice)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(nperha_rice)")

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(p2o5perha_rice)")

# Wheat sowing equation
# s(harvest_day_rice) + s(gw_2018)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(harvest_day_rice)")

# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")

# Wheat yield equation
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu4", term = "s(Res_wheat_sow)")

# Fitted values
multivariate_geo_sow_MRF_corr_endo_fitted_values <- multivariate_geo_sow_MRF_corr_endo$fitted

multivariate_geo_sow_MRF_corr_endo_fitted_values <- as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values)

# Merge the fitted results to the data and export



Irrig_Rev_rice_wheat_Mult_Res <- cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values)


summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res))

Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.574  -5.397  -0.269   4.986  44.985 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03383    3.41385   -0.01    0.992    
mu1          1.00018    0.01771   56.48   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.093 on 4525 degrees of freedom
Multiple R-squared:  0.4135,    Adjusted R-squared:  0.4134 
F-statistic:  3190 on 1 and 4525 DF,  p-value: < 2.2e-16

summary(lm(b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, Irrig_Rev_rice_wheat_Mult_Res))

Call:
lm(formula = b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, 
    data = Irrig_Rev_rice_wheat_Mult_Res)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2891 -0.7976 -0.0052  0.7454  3.3894 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.34109    0.05974   55.93   <2e-16 ***
l_ton_per_hectare  0.22766    0.01958   11.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.08 on 4525 degrees of freedom
Multiple R-squared:  0.02901,   Adjusted R-squared:  0.0288 
F-statistic: 135.2 on 1 and 4525 DF,  p-value: < 2.2e-16

plot((lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res)))

Irrig_Rev_rice_wheat_Mult_Res_sp <- Irrig_Rev_rice_wheat_Mult_Res
coordinates(Irrig_Rev_rice_wheat_Mult_Res_sp) <- c("o_largest_plot_gps_longitude", "o_largest_plot_gps_latitude")

# Map the correlations
# library(modelsummary)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- datasummary(Heading("District") * District ~ Heading("N obs") * N + Heading("%") * Percent() + lambda24 * (Mean + SD), data = Irrig_Rev_rice_wheat_Mult_Res, output = "data.frame")

# library(dplyr)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "Mean_Rice_Wheat_Rho" = "Mean")
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "SD_Rice_Wheat_Rho" = "SD")

# Irrig_Rev_rice_wheat_Mult_Res_dist <- subset(Irrig_Rev_rice_wheat_Mult_Res_dist, Irrig_Rev_rice_wheat_Mult_Res_dist$District != "Purnia")

# rice_wheat_yield_rho_dist_sf <- merge(India_aoi_sf, Irrig_Rev_rice_wheat_Mult_Res_dist, by = "District")

# rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho <- as.numeric(rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho)
# library(mapview)
# mapview(rice_wheat_yield_rho_dist_sf, zcol = "Mean_Rice_Wheat_Rho", layer.name = "Rice wheat equation correlation")
# library(sf)
# rice_wheat_yield_rho_dist_sf_sp <- as_Spatial(rice_wheat_yield_rho_dist_sf)
# library(tmap)
# tmap_mode("view")
# rice_wheat_yield_rho_dist_sf_sp_map <- tm_shape(rice_wheat_yield_rho_dist_sf_sp) +
#     tm_polygons(col = "Mean_Rice_Wheat_Rho", title = "Rice wheat equation correlation", style = "quantile") +
#     tm_layout(legend.outside = TRUE)

# tmap_save(rice_wheat_yield_rho_dist_sf_sp_map, "figures/rice_wheat_yield_rho_dist_sf_sp_map .png")

5 Factors affecting the correlation structure

f_rice_wheat_yield_MRF_corr_endo_General <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(Res_rice_sow)+s(Res_wheat_sow)+rice_duration_class_long + s(gw_2017) + s(onset_2017)+variety_type_NMWV+s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo_General <- bamlss(f_rice_wheat_yield_MRF_corr_endo_General, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)

AICc 349496.0 logPost -2543650 logLik -174245. edf 451.97 eps 1.0000 iteration   1
AICc 157756.9 logPost -2447315 logLik -78369.9 edf 457.10 eps 0.2663 iteration   2
AICc 104296.5 logPost -2420163 logLik -51655.0 edf 444.64 eps 0.2173 iteration   3
AICc 89324.62 logPost -2412631 logLik -44177.5 edf 437.79 eps 0.1257 iteration   4
AICc 86460.20 logPost -2411212 logLik -42756.4 edf 428.65 eps 0.0632 iteration   5
AICc 86239.60 logPost -2411096 logLik -42649.5 edf 425.88 eps 0.0203 iteration   6
AICc 86223.26 logPost -2411092 logLik -42646.7 edf 421.47 eps 0.0027 iteration   7
AICc 86220.64 logPost -2411088 logLik -42645.8 edf 421.15 eps 0.0008 iteration   8
AICc 86219.11 logPost -2411086 logLik -42645.3 edf 420.93 eps 0.0004 iteration   9
AICc 86218.10 logPost -2411084 logLik -42645.0 edf 420.73 eps 0.0002 iteration  10
AICc 86212.87 logPost -2411201 logLik -42644.9 edf 418.67 eps 0.0002 iteration  11
AICc 86210.00 logPost -2412487 logLik -42644.9 edf 417.53 eps 0.0002 iteration  12
AICc 86209.82 logPost -2412474 logLik -42644.9 edf 417.48 eps 0.0001 iteration  13
AICc 86209.73 logPost -2412389 logLik -42644.8 edf 417.45 eps 0.0001 iteration  14
AICc 86209.64 logPost -2412389 logLik -42644.8 edf 417.42 eps 0.0001 iteration  15
AICc 86209.64 logPost -2412389 logLik -42644.8 edf 417.42 eps 0.0001 iteration  15
elapsed time:  1.65min
Starting the sampler...

|                    |   0% 11.60min
|*                   |   5%  9.14min 28.85sec
|**                  |  10% 10.03min  1.11min
|***                 |  15%  9.29min  1.64min
|****                |  20%  9.84min  2.46min
|*****               |  25%  9.42min  3.14min
|******              |  30%  9.19min  3.94min
|*******             |  35%  8.67min  4.67min
|********            |  40%  8.11min  5.41min
|*********           |  45%  7.37min  6.03min
|**********          |  50%  6.79min  6.79min
|***********         |  55%  6.00min  7.33min
|************        |  60%  5.49min  8.23min
|*************       |  65%  4.79min  8.89min
|**************      |  70%  4.04min  9.43min
|***************     |  75%  3.42min 10.26min
|****************    |  80%  2.74min 10.97min
|*****************   |  85%  2.08min 11.81min
|******************  |  90%  1.43min 12.84min
|******************* |  95% 43.20sec 13.68min
|********************| 100%  0.00sec 14.31min

summary(multivariate_geo_sow_MRF_corr_endo_General)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo_General, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              69.2041 66.0646 69.0397 73.1000      2.030
rice_duration_class_long -1.3430 -1.9953 -1.3470 -0.6859     -1.339
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            1.283e+01 9.580e-05 5.508e-01 8.375e+01    422.419
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.897e+00 9.982e-01 1.318e+00 4.739e+00      7.444
s(onset_2017).tau21         3.733e+00 2.294e-04 1.214e-01 2.425e+01   3149.419
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.505e+00 9.995e-01 1.139e+00 3.650e+00      8.208
s(monsoon_onset_dev).tau21  1.422e+01 1.908e-03 2.439e+00 8.427e+01     15.892
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    2.514e+00 1.003e+00 2.251e+00 5.333e+00      3.656
s(median_onset_82_15).tau21 4.872e+00 2.883e-04 1.979e-02 5.178e+01     18.463
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   1.558e+00 9.998e-01 1.027e+00 4.734e+00      3.751
s(sd_onset_82_15).tau21     5.081e+00 1.665e-04 4.068e-02 4.868e+01      0.501
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       1.598e+00 9.990e-01 1.037e+00 4.329e+00      1.359
s(District,id='mrf1').tau21 1.845e+01 6.096e+00 1.642e+01 4.044e+01   1473.765
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.389e+01 3.303e+01 3.396e+01 3.444e+01     34.979
s(District,id='re2').tau21  1.636e+04 1.022e+04 1.566e+04 2.644e+04      1.087
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.599e+01 3.599e+01 3.600e+01     35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.25311 0.43654 1.14328 2.42586     -0.001
rice_duration_class_long 0.10059 0.03329 0.09880 0.16695      0.098
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             3.252e+00 1.064e-01 1.511e+00 1.579e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.642e+00 2.915e+00 5.715e+00 8.141e+00
s(g_q5305_irrig_times_rice).tau21 1.821e+00 2.364e-01 1.284e+00 6.109e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.230e+00 4.461e+00 6.231e+00 7.702e+00
s(nperha_rice).tau21              6.748e-01 5.500e-04 2.847e-01 3.440e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                3.772e+00 1.035e+00 3.737e+00 6.556e+00
s(p2o5perha_rice).tau21           1.257e-01 1.366e-04 6.628e-03 1.020e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             1.930e+00 1.010e+00 1.382e+00 4.888e+00
s(District,id='mrf1').tau21       1.431e-03 5.564e-05 7.932e-04 6.032e-03
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         1.211e+01 1.425e+00 1.140e+01 2.565e+01
s(District,id='re2').tau21        8.276e+00 2.265e+00 7.832e+00 1.717e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.584e+01 3.562e+01 3.587e+01 3.594e+01
                                  parameters
s(Res_rice_sow).tau21                 17.005
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.216
s(g_q5305_irrig_times_rice).tau21      1.765
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.581
s(nperha_rice).tau21                   1.007
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.126
s(p2o5perha_rice).tau21                0.180
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.292
s(District,id='mrf1').tau21            0.035
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             32.235
s(District,id='re2').tau21             8.974
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.893
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       110.042 107.191 110.123 112.562      3.830
variety_type_NMWV  -3.067  -3.718  -3.071  -2.438     -3.073
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     232.781    45.401   166.266   752.448   1027.121
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.951     4.433     5.909     7.452      8.387
s(District,id='mrf1').tau21   168.004    80.416   152.695   327.411   4298.841
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.867    34.741    34.870    34.961     34.994
s(District,id='re2').tau21  52603.968 33227.477 50704.077 84862.078      1.087
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.999    35.998    35.999    35.999     35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.4431753 -1.1227994 -0.4743247  0.5352021     -1.417
variety_type_NMWV    0.3437629  0.2930203  0.3434330  0.3927735      0.337
g_q5305_irrig_times  0.4227679  0.3951696  0.4233267  0.4494868      0.425
nperha               0.0015125  0.0009278  0.0015241  0.0020688      0.001
p2o5perha            0.0025713  0.0014856  0.0025638  0.0037208      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      1.273e-02 6.955e-05 1.516e-03 1.119e-01      0.179
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.455e+00 1.010e+00 1.189e+00 3.620e+00      4.475
s(District,id='mrf1').tau21 1.123e-02 3.011e-04 7.358e-03 3.616e-02      0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   2.797e+01 1.054e+01 3.035e+01 3.341e+01     31.273
s(District,id='re2').tau21  4.202e+00 9.204e-01 3.925e+00 9.290e+00      5.097
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.586e+01 3.561e+01 3.589e+01 3.595e+01     35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.216 -2.236 -2.216 -2.196     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9893 0.8978 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06595 0.04526 0.06603 0.08676      0.074
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9913 0.9271 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.146 -2.165 -2.146 -2.123     -2.141
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9892 0.9018 0.9999     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5193 0.4982 0.5194 0.5400      0.529
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9910 0.9195 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.001419 -0.018127  0.001044  0.026811     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.007514 0.003983 0.007530 0.011375      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.008788 0.003330 0.008739 0.013971      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023132 -0.008819  0.022748  0.056573      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(Res_rice_sow) + s(Res_wheat_sow) + rice_duration_class_long + 
    s(gw_2017) + s(onset_2017) + variety_type_NMWV + s(District, 
    bs = "mrf", xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                             Mean     2.5%      50%    97.5% parameters
(Intercept)              -0.17840 -0.24696 -0.17771 -0.10636      0.005
rice_duration_class_long -0.02510 -0.10378 -0.02424  0.04707     -0.019
variety_type_NMWV        -0.09329 -0.17209 -0.09272 -0.01556     -0.067
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_rice_sow).tau21       1.459e-01 1.449e-04 1.330e-02 1.129e+00      1.693
s(Res_rice_sow).alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_rice_sow).edf         1.882e+00 1.009e+00 1.487e+00 4.989e+00      5.426
s(Res_wheat_sow).tau21      1.856e-02 1.025e-04 1.995e-03 1.816e-01      0.012
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.314e+00 1.004e+00 1.083e+00 3.025e+00      1.393
s(gw_2017).tau21            5.138e-02 9.491e-05 5.993e-03 4.424e-01      0.001
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.603e+00 1.005e+00 1.267e+00 3.743e+00      1.042
s(onset_2017).tau21         2.072e-02 6.275e-05 2.181e-03 1.139e-01      1.275
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.404e+00 1.006e+00 1.183e+00 2.792e+00      4.996
s(District,id='mrf1').tau21 2.849e-03 8.893e-05 2.553e-03 7.795e-03      0.004
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.551e+01 1.640e+00 1.722e+01 2.487e+01     21.116
s(District,id='re2').tau21  7.231e-03 1.263e-04 3.328e-03 3.087e-02      0.023
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    1.101e+01 4.927e-01 9.301e+00 2.633e+01     24.188
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.06774 -0.13954 -0.07742  0.03458     -0.009
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 85816.77 logLik = -42794.54 pd = 227.6934
runtime = 862.7
---
Optimizer summary:
-
AICc = 86209.64 edf = 417.425 logLik = -42644.88
logPost = -2412389 nobs = 4527 runtime = 98.98

plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(gw_2017)")

plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(onset_2017)")

plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(Res_rice_sow)")

plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(Res_wheat_sow)")

6 Visualization

library(distreg.vis)


if (interactive()) {
   distreg.vis:: vis()
}

7 Conclusion

8 References